program pair_dist_fun;
{
CALCULATES PARTIAL PAIR DISTRIBUTION FUNCTION
INPUT  : xyz-coordinates file, specified in gdf.inp
OUTPUT : "ppdf" 
 Arbitrary (3D) unit cell allowed !
 revised Chemnitz may 97 : lots of things changed.
 Mainly - there are two different modes to run - see gdf.inp
}
const
   mxstat  = 11000;
   mxcnt   = 11000;
   dim_xyz = 20000;
   typnum  = 110;
type
   myreal = real;
   vektor = array[1..3] of myreal;
   mat3x3 = array[1..3] of vektor;
   {string = varying[50] of char;}
var
   r		    : array[1..dim_xyz] of vektor;
   drk,otp	    : array[1..dim_xyz] of integer;
   xv,xw,vk,vek,kon : vektor;
   rek		    : array[1..3] of real;
   st1,st2	    : array[0..typnum] of integer;
   nm		    : array[0..typnum] of integer;
   mynum	    : array[1..4] of real;
   mycnt	    : array[1..4] of real;
   ext		    : array[1..3] of integer;
   g,i,j,k,l,p,q,o,norm_mode,nt,dim,
   mode,aa,bb,tt,att : integer;
   V,tau,norma,alpha,cutoff,cutmin,por,acc : real;
   datapath,meno,name,mystr : string;
   txt,polet,report,error : text;
   ba:mat3x3;
   cnt:array[0..mxcnt,1..4] of integer;
   pwa:array[1..mxstat] of real;
   nnpwa:array[1..mxstat,1..4] of integer;
   npwa:integer;
   mtt:array[1..2,1..2] of integer;
   mt:array[0..typnum] of integer;
   ab:array[1..2] of integer;

   
function skalsuc(x,y : vektor):myreal;
begin
   skalsuc:=x[1]*y[1]+x[2]*y[2]+x[3]*y[3];
end;
function abshod(x:vektor):myreal;
begin
   abshod:=sqrt(skalsuc(x,x));
end;
function bound(x:myreal):myreal;
Var y:myreal;
begin
   y:=x;
   if x<-0.5 then y:=x+1;
   if x>0.5 then y:=x-1;
   if (y<-0.5) or (y>0.5) then bound:=bound(y) else bound:=y;
end;
procedure mxle3(A:mat3x3;B:vektor;Var v:vektor);
Var o,i,j:integer;
    det,dx1,dx2,dx3:real;
begin
   v[1]:=0; v[2]:=0; v[3]:=0;
   DET:=A[1,1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*A[3,1]
   +A[2,1]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*A[3,1]
   -A[1,1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*A[2,1];
   if abs(det)>0.000001 then begin
      DX1:=B[1]  *A[2,2]*A[3,3]+A[1,2]*A[2,3]*B[3]
      +B[2]  *A[3,2]*A[1,3]-A[1,3]*A[2,2]*B[3]
      -B[1]  *A[2,3]*A[3,2]-A[3,3]*A[1,2]*B[2];
      DX2:=A[1,1]*B[2]  *A[3,3]+B[1]  *A[2,3]*A[3,1]
      +A[2,1]*B[3]  *A[1,3]-A[1,3]*B[2]  *A[3,1]
      -A[1,1]*A[2,3]*B[3]  -A[3,3]*B[1]  *A[2,1];
      DX3:=A[1,1]*A[2,2]*B[3]  +A[1,2]*B[2]  *A[3,1]
      +A[2,1]*A[3,2]*B[1]  -B[1]  *A[2,2]*A[3,1]
      -A[1,1]*B[2]  *A[3,2]-B[3]  *A[1,2]*A[2,1];
      v[1]:=dx1/det; v[2]:=dx2/det; v[3]:=dx3/det;
   end;
end;
procedure readcoord;
var mrk,prk,nl:integer;
    i,j,o,p,q:integer;
begin
   for o:=0 to typnum do nm[o]:=0;
   reset(txt,meno);
   {read transformation matrix from the input file}
   for i:=1 to 3 do begin
      for j:=1 to 3 do read(txt,ba[j,i]); readln(txt);
   end;
   tt:=0;
   repeat
      tt:=tt+1;
      for q:=1 to 3 do read(txt,r[tt,q]);
      readln(txt,p);
      otp[tt]:=p;
      drk[tt]:=mt[p];
      nm[drk[tt]]:=nm[drk[tt]]+1;
   until eof(txt);
   close(txt);
   aa:=nm[1]; bb:=nm[2];
   mynum[1]:=aa+bb; mynum[2]:=aa;
   mynum[3]:=sqrt(2*aa*bb); mynum[4]:=bb;
end;
procedure writepdf_1;
var i,o	    : integer;
   x,y,r,pi : real;
   px,nx    : array[1..4] of real;
   dataf    : text;
   a	    : mat3x3;
   det	    : real;
   ii	    : integer;
begin
   for o:=1 to 4 do if mynum[o]<0.00001 then mynum[o]:=1;
   A:=ba;
   DET:=A[1,1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*A[3,1]
   +A[2,1]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*A[3,1]
   -A[1,1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*A[2,1];
   V:=abs(det);
   rewrite(dataf,'ppdf');
   o:=aa+bb;
   writeln(report,'Atoms, A, B    : ',o:1,aa:6,bb:6);
   writeln(report,'Lat_par        : ',alpha:1:3);
   writeln(report,'Atoms/volume   : ',(aa+bb):1,' ',V:1:4,' ',(aa+bb)/V:1:6);
   for i:=1 to npwa do
   begin
      write(dataf,pwa[i]*alpha:12:5);
      for o:=1 to 4 do px[o]:=1/mynum[1];
      if norm_mode<2 then 
	 for o:=1 to 4 do write(dataf,nnpwa[i,o]*px[o]:12:6);
      if norm_mode=2 then
	 for o:=1 to 4 do write(dataf,nnpwa[i,o]/2.0:12);
      writeln(dataf);
   end;
   close(dataf);
end;
procedure writepdf_n;
var i,o	    : integer;
   x,y,r,pi : real;
   px,nx    : array[1..4] of real;
   dataf    : text;
   a	    : mat3x3;
   det	    : real;
   ii	    : integer;
begin
   rewrite(dataf,'ppdfn');
   o:=aa+bb;
   for i:=1 to npwa do
   begin
      write(dataf,pwa[i]*alpha:12:5);
      for o:=1 to 4 do write(dataf,nnpwa[i,o]:12);
      writeln(dataf);
   end;
   close(dataf);
end;
procedure writepdf_0;
var i,o	     : integer;
   x,y,r,pi  : real;
   px,nx     : array[1..4] of real;
   dataf     : text;
   a	     : mat3x3;
   det,delta : real;
   ii	     : integer;
begin
   for o:=1 to 4 do if mynum[o]<0.00001 then mynum[o]:=1;
   A:=ba;
   DET:=A[1,1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*A[3,1]
   +A[2,1]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*A[3,1]
   -A[1,1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*A[2,1];
   V:=abs(det);
   pi:=3.141592654;
   delta:=cutoff/dim;
   {Note: division by dim is need for gaus.f}
   for o:=1 to 4 do nx[o]:=V/(4.0*pi*delta*delta*delta*mynum[o]*mynum[o])/dim;
   rewrite(dataf,'ppdf');
   o:=aa+bb;
   writeln(report,'Atoms, A, B    : ',o:1,aa:6,bb:6);
   writeln(report,'Lat_par        : ',alpha:1:3);
   writeln(report,'Atoms/volume   : ',(aa+bb):1,' ',V:1:4,' ',(aa+bb)/V:1:6);
   for i:=0 to dim do if cnt[i,1]>0 then
   begin
      r:=i*delta*alpha;
      x:=1.0;
      if i>0 then x:=i;
      write(dataf,r:12:4);
      if norm_mode=1 then begin
	 for o:=1 to 4 do px[o]:=nx[o]/x/x;
	 for o:=1 to 4 do write(dataf,' ',cnt[i,o]*px[o]:11:5);
      end else begin
	 for o:=1 to 4 do px[o]:=1/mynum[1];
	 for o:=1 to 4 do write(dataf,cnt[i,o]*px[o]:12:5);
      end;
      writeln(dataf);
   end;
   close(dataf);
end;
procedure pw_stat_sorted(x : real; y,pl : integer);
var i,j,o : integer;
   done	  : boolean;
begin
   if npwa=0 then begin
      npwa:=1;
      for o:=1 to 4 do nnpwa[1,o]:=0;
      nnpwa[1,1]:=pl; nnpwa[1,y]:=pl;
      pwa[1]:=x;
   end else begin
      done:=false; i:=0; repeat
	 i:=i+1;
	 if (abs(x-pwa[i])<acc) then begin
	    nnpwa[i,1]:=nnpwa[i,1]+pl;
	    nnpwa[i,y]:=nnpwa[i,y]+pl; done:=true; end;
	 if not(done) then begin
	    if (x<pwa[i]) then begin
	       j:=npwa; repeat
		  pwa[j+1]:=pwa[j];
		  for o:=1 to 4 do nnpwa[j+1,o]:=nnpwa[j,o];
		  j:=j-1;
	       until j=i-1;
	       pwa[i]:=x;
	       for o:=1 to 4 do nnpwa[i,o]:=0;
	       nnpwa[i,1]:=pl; nnpwa[i,y]:=pl;
	       npwa:=npwa+1;
	       done:=true;
	    end;
	 end;
      until (done) or (i=npwa);
      if not(done) then begin
	 npwa:=npwa+1;
	 for o:=1 to 4 do nnpwa[npwa,o]:=0;
	 nnpwa[npwa,1]:=pl; nnpwa[npwa,y]:=nnpwa[npwa,y]+pl;
	 pwa[npwa]:=x;
      end;
   end;
end; { pw_stat_sorted }
procedure compute(x,y:vektor;i,j:integer);
Var o,p1,p2,p3,mlt,bond,ii,pl:integer;
begin	   
   mlt:=1;
   for o:=1 to 3 do xw[o]:=(x[o]-y[o]);
   for p1:=-ext[1] to ext[1] do
      for p2:=-ext[2] to ext[2] do
	 for p3:=-ext[3] to ext[3] do
	    if ((i<>j) or ((i=j) and (abs(p1)+abs(p2)+abs(p3)>0))) then
	 begin
	    for o:=1 to 3 do xw[o]:=bound(xw[o]);
	    for o:=1 to 3 do
	       vek[o]:=(xw[1]+p1)*ba[o,1]+(xw[2]+p2)*ba[o,2]+
	       (xw[3]+p3)*ba[o,3];
	    por:=abshod(vek);
	    if por<acc then begin
	       writeln(error,'r = 0 for pair ',i:1,' ',j:1); end;
	    if por<cutmin then
	       writeln(error,por:1:6,i:6,j:6,otp[i]:6,otp[j]:6);
	    if por<cutoff then
	    begin
	       bond:=mtt[drk[i],drk[j]];
	       if (i=j) then pl:=1 else pl:=2;
	       if mode=1 then pw_stat_sorted(por,bond,pl);
	       ii:=round(por/cutoff*dim);
	       cnt[ii,1]:=cnt[ii,1]+pl;
	       cnt[ii,bond]:=cnt[ii,bond]+pl;
	       if npwa=mxstat then begin
		  writeln('Too many different pair distances!');
		  halt;
	       end;
	    end;
	 end;
end;
BEGIN
   mtt[1,1]:=2; mtt[1,2]:=3; mtt[2,1]:=3; mtt[2,2]:=4;
   for i:=0 to dim do
      for o:=1 to 4 do cnt[i,o]:=0;
   tau:=(1+sqrt(5))/2;
   npwa:=0;
   acc:=1e-3;
   rewrite(error,'gdf.err');
   rewrite(report,'gdf.rep');
   reset(polet,'gdf.inp');
   readln(polet,meno);
   readln(polet,alpha);
   readln(polet,cutoff);
   readln(polet,mode);
   readln(polet,cutmin);
   for o:=1 to 3 do read(polet,ext[o]); readln(polet);
   readln(polet,norm_mode,dim);
   for i:=0 to typnum do mt[i]:=0;
   readln(polet,nt);
   for i:=1 to nt do begin
      read(polet,ab[1],ab[2]);
      mt[ab[1]]:=ab[2];
      writeln('Element ',ab[1]:1,' -> ',ab[2]:1);
      if ab[2]>2 then begin
	 writeln('Projected type is <=2.'); halt;
      end;
   end; readln(polet);
   {for o:=0 to nt do read(polet,mt[o]); readln(polet);}
   close(polet);
   readcoord;
   for i:=1 to tt do if (drk[i]>0) then 
      for j:=i to tt do if (drk[j]>0) then compute(r[i],r[j],i,j);
   if mode=1 then begin writepdf_1; writepdf_n; end else writepdf_0;
   close(report);
   close(error);
END.

